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17 Abstract 

18 A number of social-ecological systems exhibit complex behaviour associated with nonlinearities, 

19 bifurcations, and interaction with stochastic drivers. These systems are often prone to abrupt 

20 and unexpected instabilities and state shifts that emerge as a discontinuous response to gradual 

21 changes in environmental drivers. Predicting such behaviours is crucial to the prevention of or 

22 preparation for unwanted regime shifts. Recent research in ecology has investigated early 

23 warning signs that anticipate the divergence of univariate ecosystem dynamics from a stable 

24 attractor. To date, leading indicators of instability in systems with multiple interacting 

25 components have remained poorly investigated. This is a major limitation in the understanding 

26 of the dynamics of complex social-ecological networks. Here, we develop a theoretical framework 

27 to demonstrate that rising variance - measured, for example, by the maximum element of the 

28 covariance matrix of the network - is an effective leading indicator of network instability. We 

29 show that its reliability and robustness depend more on the sign of the interactions within the 

30 network than the network structure or noise intensity. Mutualistic, scale free and small world 

31 networks are less stable than their antagonistic or random counterparts but their instability is 

32 more reliably predicted by this leading indicator. These results provide new advances in 

33 multidimensional early warning analysis and offer a framework to evaluate the resilience of 

34 social-ecological networks. 
35 

36 
37 
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38 Introduction 

39 Social-ecological systems are often difficult to investigate and manage because of their inherent 

40 complexity (1). Small variations in external drivers can lead to abrupt changes associated with 

41 instabilities and bifurcations in the underlying dynamics (2-4). These transitions can occur in a 

42 variety of ecological and social systems, and are often unexpected and difficult to revert (4). 

43 Anticipating critical transitions and divergence from the present state of the system is 

44 particularly crucial to the prevention or mitigation of the effects of unwanted and irreversible 

45 changes (5-10). Recent research in ecology has focused on leading indicators of regime shift in 

46 ecosystems characterized by one state variable (5,7,11,12). These indicators are typically 

47 associated with the critical slowing down phenomenon: as the system approaches a critical 

48 transition, its response to small perturbations of the stable state becomes slower (11). It has 

49 been shown that in univariate systems (i.e., with only one state variable) critical slowing down 

50 entails an increase in the temporal variance and autocorrelation of the state variable (5). The 

51 case of systems with several mutually interacting components, however, has remained poorly 

52 investigated (13-15), while the connection between network stability and research on indicators 

53 for loss of resilience has been elusive (16). 

54 Here we develop a theoretical framework to analyze early warning signs of instability and regime 

55 shift in complex networks. We provide analytical expressions for a set of precursors of instability 

56 in complex systems with additive noise for a variety of network structures. 

57 We consider a social-ecological system with N components (nodes) coupled through a set of 

58 links. The state of the system is expressed by the vector x of length N, whose terms x, represent 

59 the state of node i. The local stability of a state x* is evaluated through a linearization, — = A • y, 

60 where y=x-x* is the displacement of x from x*; A is the JVxJV matrix expressing the interactions 
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61 among nodes in the (linearized) dynamics (see Methods). In population ecology this framework 

62 is typically used to express the dynamics of a community of N populations interacting according 

63 to the relationships determined by the matrix A, often known as "community matrix" (2,17 -19); 

64 likewise, in social systems A describes the network of interactions (e.g., trade, migration, flow of 

65 information among people, groups of individuals, or countries (20-23)). The off-diagonal terms of 

66 A determine the pairs of interacting nodes as well as the strength of their interaction. The 

67 dynamics are stable if the maximum real part of the eigenvalues of A, Max[Re{X}], is negative. 
68 

69 Classic ecological theories (2,3) have considered the case of networks with randomly connected 

70 nodes (with a certain probability, Q. The strength (p) of the interactions between them is 

71 represented by a zero-mean random variable of variance a 2 . May (2,3) showed that random 

72 networks become unstable as connectivity (i.e., Q, size (i.e., M) or the strength variance increase. 

73 The stability of networks with prescribed architectures (e.g., predator-prey, competitive or 

74 mutualistic interactions) also depends on connectivity, strength variance, system size, as well as 

75 on the network structure (17,19). 

76 More generally, the off-diagonal terms of A may result from a set of "rules" expressed as a 

77 function of a few parameters of which connectivity and strength variance are just an example. 

78 Changes in the structure and intensity of the interactions correspond to variations in these 

79 parameters, which, in turn, can lead to instability by modifying the community matrix and its 

80 eigenvalues. How can we evaluate whether ongoing changes in the interactions within a social- 

81 ecological network are reducing its resilience? Is there a way to use measurable quantities to 

82 determine whether the system is about to become unstable? 
83 
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84 In one-dimensional systems leading indicators are typically associated with behaviors resulting 

85 from the eigenvalue tending to zero at the onset of instability. This effect entails a slower return 

86 to equilibrium after a "small" perturbation (11,24). Known as "critical slowing down", this 

87 phenomenon exists also in systems with multiple interacting components, though it is hard to 

88 recognize and therefore it does not constitute an effective leading indicator of instability. In fact, 

89 in "real world" applications the equations driving the dynamics are not known and, therefore, the 

90 network nodes in which slowing down is expected to occur are not known a priori. Critical 

91 slowing down, however, has been related to an increase in variance and autocorrelation in the 

92 state variable of one dimensional systems (5,7,25). Here we provide a theoretical framework to 

93 investigate early-warnings in the variance, autocorrelation, and power spectrum of multi- 

94 dimensional systems with interactions described by a given network structure. 
95 

96 Methods 

97 We consider a network with N interacting nodes. The state of the system, x={xi, xz, ... xn} , is 

98 governed by dynamics: dx = f(x,p,C~)dt + v I dW , where f-{fi, fi, fn} is a iV-dimensional 

99 vector function expressing the deterministic component of the dynamics of x, as a function of a 

100 set of parameters, p and C; I is the identity matrix, and v dW is an additive stochastic driver 

101 represented by a white Gaussian noise of mean zero and intensity vdt. If we consider a small 

102 perturbation y forcing the system away from its equilibrium point x* (i.e., y=x-x*), inserting x= 

103 x*+y in the above equation and linearizing f[x*+y, p, Q around x* we obtain 

104 dy = A(j?)ydt+IvdW , (l) 

105 

106 where A. = \d f i I dXj\ =J * . Eq. (1) is a multivariate Ornstein-Uhlenbeck process (26). 
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107 The stable states, x*, of Eq. (1) are the same as those of their deterministic counterparts, 

108 ^ = f(x,p,C~) (27). These states are stable if the maximum real part of the eigenvalues of A is 

109 negative. To identify early warning signs of network instability, we relate the steady state 

110 covariance matrix S y = (y s , y T s ) to the eigenvalues, X, of A, where y s is calculated from the steady 

111 state solution of Eq. (1). We first look for leading indicators of instability in the behavior of the 

112 covariance matrix, S y , of y as the system approaches instability. The (/,/) element of S y is: 

113 S y (i,j) = {y j y j ^ — {y i )^y j ^, where ( ) represents the average. The covariance matrix of the 

114 stationary dynamics of the system can be obtained (26) as the solution of Eq. (2): 
115 

116 A(p, C)S y + S y A T (p, C) = -v/ (2) 
117 

118 S y is a function of the linearization matrix, A[p,C), which, in turn, depends on the control 

119 parameters [p or C). At the onset of instability (i.e., as Max[Re(X)] — >0) the maximum element of 

120 the covariance matrix, S y , of y increases. More details on the time-lag correlation and power 

121 spectrum can be found in the Supplementary Materials. While the linearization matrix, A, here 

122 accounts the interconnections existing among nodes within the network (i.e. the pairs of nodes 

123 that are connected by a link (3, 17)), the covariance matrix, S y , expresses the variance of the 

124 fluctuations of the state variable at each node (diagonal terms) and the interrelationship 

125 (positive or negative) of the fluctuations between pairs of nodes (off-diagonal terms). To better 

126 understand the structure of S y , we look at the case of a network with only two nodes. In this case 

127 the above equation for the covariance matrix can be solved analytically, and the covariance 

128 matrix reads (26) 
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0 , . -vZ)e?(A)I + (rr(A)I-A)vI(7V(A)I-A) r 

129 S(p) = , 

} F 2Det(A)Tr(A) 

130 where Det is the determinant and Tr the trace of the matrix, which can be expressed as a function 

131 of the eigenvalues, Ai,2. 

132 Thus, the covariance matrix diverges for Xi, A.2— >0. The time correlation matrix, p y (A), can also be 

133 computed analytically and also diverges for X\, ta— >0, independently of A. The general analytical 

134 expressions of p y (A) and of the power spectrum of y are also reported in the Supplementary 

135 Materials. 

136 We generate networks of size N, with a variety of architectures for A (see Supplementary 

137 Materials), and reach instability either by keeping constant the connectivity, C, while changing 

138 the strength of the interactions, p, or by varying C for a fixed p (2,17,19). We then use the 

139 analytical relationship between the steady state covariance matrix, S y , of y and the eigenvalues of 

140 the matrix A (Eq. (2)). Similarly, we express the time-lag correlation, p y , and the power spectrum, 

141 P y , of y as a function of A and its eigenvalues. 
142 

143 Results and discussion 

144 We find that the elements of both S y and p y increase as the system approaches instability (i.e., 

145 Max[Re(A)]— >0). Therefore, we investigate potential indicators for early warning in the behavior 

146 of suitable components of S y , p y and P y for Max[Re(A)]— >0. To that end we first consider the 

147 components of S y corresponding to the most connected, the most central (28) and the least 

148 connected nodes of the network. We also consider indicators based on the properties of the 

149 entire network, such as the maximum and the difference between the maximum and minimum of 

150 the matrix S y . 

7 



Downloaded from http://biorxiv.org/on September 18, 2014 

151 Most of the indicators based on the covariance matrix, S y , have a non-trivial dependence on 

152 Max[Re(X)] (see Figures 1, S1-S6). The maximum element of S y (Max[S y ]) and Max[S y ]-Min[S y ] 

153 provide the most effective indicator of early warning in most networks (Figures 1, S7 and S8). In 

154 mutualistic (++) networks Max[S y ] corresponds to the most connected node (the "hub"), 

155 regardless of their topological structure (Supplementary Materials, Figures S9-S10). All these 

156 indicators based on S y improve their performances when the size, N, of the network increases 

157 (compare main panels to insets in Figure 1; see also Supplementary Materials, Figure Sll). Thus 

158 our ability to detect early warning signs and predict tipping points is enhanced in more diverse 

159 systems (16). 

160 We also look at the relationship between the maximum element of the time-lag correlation 

161 matrix, p y (A) (where A is the time lag), and Max[Re(X)] for different values of A, p and C (Figures 

162 S12-S14). Although significant, these indicators are less efficient with respect to the case with 

163 zero time-lag (i.e., indicators based on S y ). Finally, the power spectrum does not appear to be an 

164 effective indicator, as we identified only weak changes in P y for increasing values of p and 

165 Max[Re(X)] (see Supplementary Materials, Figure S15). Therefore, here we focus on early 

166 warning signs provided by the way Max[S y ] varies as a function of changes in Max[Re(X)] . 

167 A warning sign is effective if (a) it appears in time to prevent (or prepare for) the occurrence of 

168 instability (29-30); (b) it relies on a well-defined and easy to recognize indicator (e.g., a 

169 detectable or significant increase in variance (29-30)); and (c) it does not give false positives (or 

170 false negatives) (31). We use these criteria to evaluate the effectiveness of Mc?x[S y ] as a leading 

171 indicator of instability with different network structures and levels of noise (32). 

172 To investigate the effect of noise, we first consider the "mean-field" case of networks in which the 

173 absolute value of the interaction strength between connected nodes is a constant, p; we gradually 
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174 increase p or C until Max[Re(X)] becomes positive (17-19). We observe (Figure 2) a consistent 

175 increase in Max[S y ] for all network structures, regardless of whether instability is attained by 

176 increasing interaction strength or connectivity (Figures S1-S6). The network structure, however, 

177 affects the timeliness of Max[S y ] as a leading indicator. In fact, Max[S y ] exhibits a more defined 

178 increase and a better anticipation of the onset of instability in the case of random networks than 

179 with all the other structures. In the case of these "mean field" networks we did not consider the 

180 antagonistic structure because antagonistic networks with constant interaction strength (in 

181 absolute value) are always stable regardless of the parameters p and C (see Supplementary 

182 Materials, Figs. SI, S2). 

183 Likewise, in the case of random interaction strengths Max[S y ] exhibits a well-defined increase 

184 and a better anticipation of the instability in random networks than with the more organized 

185 structures typical of ecological or social systems (Figures 3, S3-S6). The seemingly weaker 

186 increase in Max[S y ] observed in the social-ecological networks is only an apparent effect of the 

187 scale. Indeed, as it will be shown later, suitable detection criteria of early warnings are more 

188 successful in mutualistic networks than in their random counterparts. Noise has the effect of 

189 amplifying the intensity of the warning sign (compare the scales in Figs. 2 and 3), while inducing 

190 weak random fluctuations with no substantial impact on the overall behavior of Max[S y ] at the 

191 onset of instability (see Supplementary Materials). In scale free networks the increase in Max[S y ] 

192 (Figure 3) is again only apparently muted. In fact, in these networks detection criteria are quite 

193 successful in recognizing early warning signs (Figure 4); moreover, local indicators (e.g., the 

194 variance of the most central node) can exhibit a more pronounced increase that can be used as 

195 an early warning sign of instability (Figure 1 and Supplementary Materials). 
196 
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197 Conclusions 

198 We have identified some suitable early warning signs in social-ecological networks in agreement 

199 with those identified by Ref. (33), and provided a theoretical framework for their interpretation. 

200 Overall, the performances of Max[S y ] as a leading indicator of instability change between random, 

201 antagonistic, mutualistic/social networks. This indicator gives an earlier and "sharper" warning 

202 sign in random than mutualistic and social networks. The warning sign, however, is harder to 

203 detect and is more likely to be missed in random and antagonistic networks than in their 

204 mutualistic or social counterparts (Figures 4, S16). Thus, by affecting the probability that early 

205 warnings are missed, the sign of the interactions within the network determines the consistency 

206 and reliability of this leading indicator. In fact, different realizations of the same network 

207 dynamics can yield different results in the behavior of Max[S y ] and thus this indicator might not 

208 detect in useful advance the emergence of instability (Figures S17-S18). The probability of true 

209 positives is close to 100% (i.e., negligible probability of false negatives) in mutualistic networks, 

210 and much smaller in random and antagonistic (predator-prey, cascade or compartment) 

211 networks (Figures 4, S16). Thus, while mutualistic networks are less stable than their 

212 antagonistic counterparts (17), their instability can be predicted with less uncertainty. An 

213 increase in Max[S y ], however, would not provide information on how close the system is to the 

214 onset of instability. Rather, it would just indicate that the system is losing resilience and 

215 approaching unstable conditions (29). Therefore, in contrast to previous expectations (16), it is 

216 not the heterogeneity in the topology of the network that plays a key role in the abruptness of 

217 critical transitions and our ability to predict them. Rather, it is the type of interactions between 

218 the nodes that determines how networks respond to external perturbations. In fact, there is a 

219 trade-off between local and systemic resilience: mutualism (++) is associated with a reduced 
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220 local stability and resilience of the system (17,19), but does not induce abrupt critical transitions. 

221 In contrast, networks with mixtures of interaction types (+-,++, --) exhibit shorter recovery times 

222 after displacement from equilibrium (i.e., a stronger local resilience) (17-18), but in these 

223 systems the emergence of systemic instability and critical transitions is more difficult to predict 

224 in useful advance. 

225 This study combines stability theories from community ecology (2,17) to recent research on 

226 indicators of critical transition (7,9,16), and develops a unified framework that offers a new 

227 perspective for the evaluation of the resilience and anticipation of instability in social-ecological 

228 networks. 
229 
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310 Figure 1. Leading indicators of instability based on different elements of the covariance matrix 

311 (S y ), including the maximum (in absolute value) element, Max[S y ], the difference between 

312 Max[S y ] and Min[S y ], the element of S y corresponding to the most connected, least connected, or 

313 highest eigenvector centrality (24) network node. Random [left) and scale free (right) (30) 

314 network generated with N-SO and C=0.1 (main panels) and N-0.1 and C-0.S (insets). Instability 

315 (i.e., decrease in Max[Re[X)]) is attained by increasing the interaction strength p (mean field 

316 case). The figures represent average behavior over 100 realizations. 

317 
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321 Figure 2. Max[S y ]as a leading indicator of instability in a "mean field" network with constant 

322 interaction intensity (in absolute value) , p. Instability is attained by increasing p (main panel A, 

323 with N-20, C-0.2) or C (inset B, with N-20, and C increasing from 0.1 to 1) with different network 

324 structures. The figures represent average behavior over 1000 realizations. 

325 
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332 Figure 3. A) Case with random interaction strength (see methods). Main panel: instability is 

333 reached by increasing p (with N=20; C-0.2). First inset (B): p is constant while C increases 

334 between 0.1 and 1. C) Same as the first inset (B) but only for the scale-free network (notice the 

335 different scale on the vertical axis). The figures represent average behavior over 1000 

336 realizations. 
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338 

339 Figure 4. Distribution of the correlation, pK, between Max[S y ] and the parameter p, after 1000 

340 realizations for the full disordered (not mean-field) case. If pk is significant (p-value<0.05) and pk 

341 > 0.5 the increase in Max[S y ] is interpreted as an early warning sign. We calculate these detection 

342 statistics for several realizations of each network structure and determine the probability of 

343 detecting the early warning sign of instability. We consider eleven different network 

344 architectures typical of ecological or social networks, including random (R), predator-prey (PP), 

345 cascade (Case), compartmentalized (Comp), mutualistic (M), bipartite (Bip), nested (N), nested 

346 with competition (N+C), scale free (SF), and small world (SW). These networks have different 

347 structures for the adjacency matrix and different combination of interaction types, i.e (++) 

348 mutualistic, (+-) antagonistic, (--) competitive or a combination of them (See Supplementary 

349 Materials for more details). 

19 
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1 Supplementary online materials 

2 SI. Structures of Socio-Ecological Interacting systems. 

3 a) Random matrix (1) with connectivity C. We pick an element Ay at random, and with 

4 probability C we assign a link between nodes i-j and ;'-/' (Ay). Each of these two links has a 

5 weight (or strength), p, that is positive (or negative) with probability 0.5; b) Antagonistic 

6 matrix (2), where connections are again random (with probability Q, but if Ay has sign +(-), 

7 then Ay has sign - (+); c) Cascade network (3): links occur with a given probability, but form a 

8 hierarchical structure, whereby there is a top "predator" (with sign +) that feeds on all the 

9 other species (with sign -); then there is the second top predator, and so on; species left with 

10 the lowest ranking function as producers; d) Compartmentalized structure (5), formed by 

11 groups of antagonistic species that interact only with species within their own group; e) 

12 Mutualistic network (6), with a random matrix, but where both Ay and Ay have positive signs 

13 (++); f) Modular mutualistic matrix (5): nodes are divided into communities that positively 

14 interact within their groups (++). g) Bipartite Mutualistic matrix (6): nodes are divided into 

15 two groups, and each node positively interacts only with nodes of the other group (++); h) 

16 Bipartite Nested Mutualistic (7): a bipartite graph is generated with hierarchical structure 

17 where specialist nodes (i.e., with only few mutualistic links) tend to interact with a suitable 

18 subset of the mutualistic partners of the generalist nodes (++); i) Nested Mutualistic with 

19 Competition (8): nodes are divided into two groups; each node positively interact with nodes 

20 of the opposite group (++), while compete with nodes of the same group (--); j) Barabasi- 

21 Albert (BA) networks (9): binary one-zero networks are generated with power-law degree 

22 distribution. This algorithm simulates a preferential attachment process, in which a new 

23 vertex with d edges is added at each step. The BA graph displays a scale-free behavior that 

24 strongly correlates with the network's robustness to failure. We then assign to each link a 

1 



Downloaded from http://biorxiv.org/on September 18, 2014 

25 positive weight p; The connectivity of the BA networks is controlled indirectly by the 

26 parameter d. In our simulation we set d-1 to generate low connectivity networks, d-2 for 

27 average connectivity and d-A for high connectivity, k) Watts-Strogatz (WS) network (10): this 

28 is a one-parameter (degree of disorder r) model that interpolates between an ordered finite 

29 dimensional lattice and a random graph. Its main property is to display both high clustering 

30 coefficient and small world property. The WS model displays this duality for a wide range of 

31 the rewiring probabilities r. In our simulation we used r=0.3. The connectivity of the WS 

32 model depends on the starting 2k-regular graph (10). For low connectivity we set k=2; for 

33 average connectivity k=4; for high connectivity k>5. 

34 The linkage density in each network is described by the parameter C (the connectivity ). 

35 We generate each of the structures (a)-(/c) keeping constant the connectivity C and changing 

36 the strength of the interactions (alternatively, we fix p and vary Q, until the dynamics become 

37 unstable (1,6). Thus, we use the parameter, p, as an indicator of the weights of the interactions 

38 between connected nodes (connectivity). We consider three different levels of disorder in the 

39 control parameter p: 1) Mean Field: where we assign a constant value, p. We do not consider 

40 antagonistic mean field networks (i.e., with constant |p| but (+-) interactions) because they 

41 are always stable (see below); 2) Weak disorder, the interaction strength is drawn from a 

42 normal distribution with mean p and standard deviation 0.1 p (i.e., JV(p, 0.1 p); 3) Strong 

43 disorder, interaction strengths are drawn at random from ±\N [0, p)\ (6). Therefore, in all 

44 these cases, by increasing the value of p (C) the resilience decreases until the system become 

45 unstable for p= p c , or C=C C (6). 

46 
47 
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48 S2.1 Early Warnings in Multidimensional Mean Field Networks 

49 

50 In the mean field case with (++) interactions, the matrix is symmetric (A=A T ) and therefore all 

51 eigenvalues are real numbers. In this case we note that, if U is the matrix of the eigenvectors 

52 of A (i.e. A Ui= Aj Ui), then U A U T =diag(A), where diag( X) is a diagonal matrix with all the 

53 eigenvalues of A. Therefore in this case we can write Eq. (2) in the main text as: 

54 UAU T US y U T +US y U T UAU T =-vI =diag{AyL + Ydiag(A) , 

55 where UU T = I and £ = U S y U T is a diagonal matrix whose z'-th eigenvalue is Q = — — . 

2 X i 

56 Therefore if the socio-ecological system reaches an instability (i.e. h=0) for a given /, then 

57 there is at least one element of the matrix £ and thus of S y that diverges. 

58 We also notice that in the mean field case we omit the analysis of early-warning performed on 

59 antagonistic (+-) matrices, A'. In fact, if A' is anti-symmetric (i.e., A'=-A' T , and thus with 

60 predator prey, cascade or compartmentalized interactions), then the early warning analysis 

61 on those matrices is trivial as their stability does not depend on p or C. In fact, if we multiply 

62 the off-diagonal terms of those matrices by the imaginary unit /, i.e. B-i (A'-I A'), then B is 

63 Hermitian and has all real eigenvalues. Therefore A'-I A' has all pure complex eigenvalues 

64 independently of p and C, which means that A' is always stable - given that the self interaction 

65 terms are negative constant, i.e. A'i,i=-d<0. 

66 Figures SI and S2 show a rise in Max[S y ] as the system approaches instability in mean field 

67 networks. 
68 

69 
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Figure SI. Increase in Max[S y ] as Max[Re(X)] tends to zero for mean field networks of size, 
N-20, C=0.2. Increasing values of Max[Re(X)] are obtained by increasing the interaction 
strength, p. The plotted values are the ensemble averages of 1000 realizations. 
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Figure S2. Increase in Max[S y ] as Max[Re(A)]->0 for mean field networks of size, N-20, p«p c . 
Increasing values of Max[Re(X)] are obtained by increasing the connectivity, C. The plotted 
values are the ensemble averages of 1000 realizations. 
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85 S2.2 Early Warnings in Complex Networks with Weak Disorder 
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Figure S3. Increase in Max[S y ] as Max[Re(X)] tends to zoro for complex networks with "weak" 
disorder (see Section 1) of size, N-20 and C=0.2. Increasing values of Max[Re(X)] are obtained 
by increasing the interaction strength, p. The plotted values are the ensemble averages of 
1000 realizations. 
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97 Figure S4. Increase in Max[S y ] as Max[Re(A)]->0 for complex networks with "weak" disorder 

98 (see Section 1) of size, N-20 and p«p c .. Increasing values of Max[Re(X)] are obtained by 

99 increasing the connectivity, C. The plotted values are the ensemble averages of 1000 
100 realizations. 
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S2.3 Early Warnings in Complex Networks with Strong Disorder 
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Figure S5. Increase in Mc/x[S y ] as Max[Re(X)] tends to zero for complex networks with 
"strong" disorder (see Section 1) of size, N-20 and C-0.2. Increasing values of Max[Re(X}] are 
obtained by increasing the interaction strength, p. The plotted values are the ensemble 
averages of 1000 realizations. 
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Figure S6. Increase in Max[S y ] as Max[Re{X)] tends to zero for complex networks with 
"strong" disorder (see Section 1) of size, N-20. Increasing values of Max[Re(X)] are obtained 
by increasing the connectivity, C. The plotted values are the ensemble averages of 1000 
realizations. 

S3. Indicators and Node Properties 

In this section we show more analysis on the effectiveness of different node variances as 
precursors of instability. 
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Figure S7. Elements of the covariance matrix S y corresponding to nodes with the highest 
number of connections (green), lowest number of connections (light blue), highest 
eigenvector centrality (gold), max[S y ] (violet) and max[S y ]-min[S y ] (purple), in the case of: 
(A) mutualistic, (B) mutualistic nested with competition, (C) small world interactions, for 
mean field networks (of size N-20 and connectivity C-0.3). The plotted values are the 
ensemble averages of 1000 realizations. 
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Figure S8. Elements of the covariance matrix S y corresponding to nodes with the highest 
number of connections (green), lowest number of connections (light blue), highest 
eigenvector centrality (gold), max[S y ] (violet) and max[S y ]-min[S y ] (purple) in the case of (A) 
random, (B) predator-prey, (C) mutualistic, (D) mutualistic nested with competition, (E) Small 
world, (F) Barabasi-Albert, networks with strong disorder (of size N-20 and connectivity 
C=0.3). The plotted values are the ensemble averages of 1000 realizations. 
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Figure S9. Frequency distribution of the degrees (i.e., number of connections) of the 
networks' nodes and (with partially filled circles) of the nodes associated with the maximum 
value of the covariance matrix S y in mean field networks with a variety of interactions. Based 
on a set of 100 realizations. Notice how, in mutualistic networks the node corresponding to 
max[S y ] is associated with the nodes with the highest degrees (i.e. the generalist species). 
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149 

150 Figure S10. Frequency distribution of the degrees (i.e., number of connections) of the 

151 networks' nodes, and (with partially filled circles) frequency distribution of the nodes 

152 associated with the maximum value of the covariance matrix S y in "strongly" disorganized 

153 networks with a variety of interactions. Based on a set of 100 realizations. Notice how, in 

154 mutualistic networks the node corresponding to max[S y ] is associated with the nodes with the 

155 highest degrees (i.e. the generalist species). 
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159 
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1,00 




161 

162 Figure Sll. Effect of the network size on the magnitude of the early warning sign. Maximum 

163 real part (in absolute value) of the network's eigenvalues as a function of the network size, N 

164 for a random network with C=0.25 and p=p c =l/VNC. As N increases the max of Re(A) tends to 

165 zero as Max[Re(X)] ~ N- L5 and the resilience of the system decreases, while the "height" of the 

166 early warning increases. Therefore, as N increases, the early warning sign becomes sharper 

167 (see also Figure 1 in the main text). 
168 
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S4. Early Warning for Time Correlation Matrix and Power Spectrum 



We calculate the time lag correlation matrix, p y (A) (where A is the time lag), and the power 
spectrum, P y (co) (where a) is the frequency) for the steady state dynamics: 



p(A) = (y f (f).y/ (t + A) = a(p)Exp{A{p)A) 



(3) 



(iQ)-A(p)y l vI(i(o + A I (p)) 



In 



(4) 



The following figures show the behavior of the autocorrelation and power spectrum as the 
system approaches instability. 
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Figure S12. Increase in Max[p y ] as Max[Re(X)] tends to zero for "strongly" disordered 
networks with a random architecture with N-20 and C-0.3. Increasing values of Max[Re[X)] 
are obtained by increasing the interaction strength, p. The plotted values are ensemble 
averages of 100 realizations. 
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186 Figure S13. Increase in Max[p y ] as Max[Re(X)] tends to zero for "strongly" disordered 

187 networks with a predator-prey architecture and N-20, C-0.3. Increasing values of Max[Re(X)] 

188 are obtained by increasing the interaction strength, p. The plotted values are ensemble 

189 averages of 100 realizations. 
190 
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Figure S14. Increase in Max[p y ] as Max[Re(X)] tends to zero for "strongly" disordered 
networks with a mutualistic architecture and N-20, C-0.3. Increasing values of Max[Re(X}] 
are obtained by increasing the interaction strength, p. The plotted values are ensemble 
averages of 100 realizations. 
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Figure SI 5. Left panel: max element of the power spectrum matrix as a function of frequency 
for three different architectures. The impact of the structure on the spectrum is negligible. 
Right panel: Power spectrum evaluated in the minimum and maximum frequency as p tends 
to p c (and thus Max[Re(X)] tends to zero) for strongly disordered systems (N-20 and O0.2) 
with random, predator-prey, and mutualistic interactions. Increasing values of Max[Re(X}] 
lead to a decrease in 5= Max[P(co m i n )] - Max[P(co max )], that therefore might be considered a 

206 precursor for a critical transition. However, the intensity of this early warning sign is quite 

207 weak and thus difficult to detect The plotted values are ensemble averages of 100 

208 realizations. 
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217 S5. Early Warning Detection 

218 To evaluate whether the onset to instability can be anticipated in time by an increase in 

219 Max[S y ] (or in other suitably chosen elements of S y ), we test the correlation (11) between 

220 Max[S y ] and the control parameter (p or Q that is gradually varied to increase Max[Re(X)] up 

221 to a given threshold (here chosen equal to -0.2). If the correlation, pk, (evaluated with the 

222 Kendall-x test) is significant and greater than 0.5, the increase in Max[S y ] is interpreted as an 

223 early warning sign. We repeat this analysis for 1000 realizations of the random interaction 

224 strength network and determine the distribution of correlations along with the number of 

225 realizations with positive warning sign. 

226 
227 
228 
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230 Figure S16. Probability of detecting true positives (i.e. of not missing a warning sign) in the 

231 case of mean field networks, using the same detection criteria as in Figure 4. 
232 
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238 

239 
240 
241 



18 



Downloaded from http://biorxiv.org/on September 18, 2014 



242 



Random 



Q 
Q_ 



m 






rrr>^ 



Predator— Prey 



Cascade 



0.50 0.55 0.60 0.65 

Kendalh 



Compartimentalized 




15 

LL 
CL 

5 
0 



0.55 0.60 0.65 

Kendalh 

Mutualisitc 



0.50 0.52 0.54 0.56 0.58 0.60 0.62 

Kendalh 

Modular 




Q 
Q- 



20 
15 

LL 

Q 10 
CL 



0.50 0.52 0.54 0.56 0.58 0.60 



Kendalh 

Bipartite 




1 

om 

0 .503 .550 .60) .650 .700 .750 .800.85 

Kendalh 

Nested 





0.85 0.90 0.95 

Kendalh 

Nested with Competition 



0.80 0.85 0.90 0.95 

Kendalh 



Scale Free 



Q 4 
CL 

2 
0 



0.5O0.55O.6OO.65O.7OO.75O.8O0.85 

Kendalh 

Small World 



0.50 0.55 0.60 0.65 0.70 

Kendalh 





0.88 0.90 0.92 0.94 0.96 

Kendalh 



0.88 0.90 0.92 0.94 0.96 0.98 

Kendalh 



243 Figure S18. Frequency distribution of the pK statictics used to detect early warning signs of 

244 instability in the case of full disordered networks. 
245 
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